nothing 



Energy-Dependent Harmonic Ratios of the Cyclotron Features of 
X0331+53 in the 2004-2005 Outburst 

M. Nakajima^'^, T. Mihara^ and K. Makishima^'^ 
ABSTRACT 

We report on changes of the cyclotron resonance energies of the recurrent 
transient pulsar, X033H-53 (V0332+53). The whole RXTE data acquired in 
the 2004-2005 outburst were utilized. The 3 — 80 keV source luminosity varied 
between 1.7 x 10^^ and 3.5 x 10^* erg s~^, assuming a distance of 7 kpc. We 
confirmed that the fundamental cyclotron resonance energy changed from ~ 22 
to ~ 27 keV in a clear anti-correlation to the source luminosity, and without 
any hysteresis effects between the rising and declining phases of the outburst. In 
contrast, the second harmonic energy changed from ~ 49 to ~ 54 keV, implying a 
weaker fractional change as a function of the luminosity. As a result, the observed 
resonance energy ratio between the second harmonic and the fundamental was 
~ 2.2 when the source was most luminous, whereas the ratio decreased to the 
nominal value of 2.0 at the least luminous state. Although the significance of 
this effect is model dependent, these results suggest that the fundamental and 
second harmonic resonances represent different heights in the accretion column, 
depending on the mass accretion rate. 

Subject headings: pulsars: individual(X033H-53, 4U0115-I-63) — X-rays: bina- 
ries 
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INTRODUCTION 



Magnetic field strengths on the surface of accreting X-ray pulsars can be measured accu- 
rately by observing Cyclotron Resonant Scattering Feature (CRSF), because the fundamental 
electron cyclotron resonance energy Ea_i is described as Eai = 11.6 Bi2{l + Zg)~^ (keV). Here, 
Bi2 is the magnetic field strength in units of 10^^ Gauss, and Zg is the gravitational redshift. 
Applying this relation to X-ray detected CRSFs, the magnetic field strengths of ~15 X-ray 



1995: Makishima et al. 1999; Coburn et al. 



pulsars have been accurately measured (iTriimper et al. I Il978l : IWhite et al. I Il983l : iMihara 



2002 



and re ferences therein). The r esults are 



clustered in a relatively narrow range of (1 — 4) x 10^^ G (IMakishima et al. Ill999l ). 



Since the cyclotron resonant scattering occurs near the pulsar surface, the resonance 
energies were believed to be constant and intrinsic to each pulsar. However, an unexpected 
CRSF energy change was found from the recurr ent transient pulsar 4U 0115+63 w i th th e 
Ginga observations performed in 1990 and 1991 ( Mihara Ill995l : iMihara et al. Ill998l 120041 ): 
in the 1991 minor outburst when the luminosity was ~ y of that in the 1990 and other typical 
outbursts of this object, the value of £^ai was ~ 1.4 times higher than the so far reported 
i?ai~ 11 keV. In order to interpret the change in Ea,i, Mihara et al. (2004) proposed that the 
height of cyclotron scattering region in the accretion column changes depending on the X-ray 
luminosity. A further study of the luminosity- dependent CRSFs change was conducted by 
Nakajima (2006), Nakajima et al. (2006a,b), and Tsygankov et al. (2007), using Rossi X-ray 
Timing Explorer {RXTE) data of 4U 0115+63 acquired in 1999 and 2004. According to 
these results, E^i increased from ~ 10 to ~ 16 keV as the source luminosity decreased across 
a narrow range of (2 — 4) x 10^^ erg s~^ (at an assumed distance of 7 kpc; Negueruela & 



Okazaki 2001). In addition, the second harmonics observed at 
started to change. 



20 keV disappeared as i?ai 



This luminosity-dependent change in the CRSF energy has been found from another 
source, X0331+53 (V0332+53), with Ginga, INTEGRAL and i^ZTE observations (Mihara 
et al. 1998; Mowlavi et al. 2006; Nakajima 2006; Tsygankov et al. 2006, 2009). This source 
is a recurrent transient pulsar, and exhibits a very prominent fundamental CRSF at 



28 keV of the X-ray spectrum, as first discovered with G inga (IMakishima et al. I Il990al ) 



after a suggestion with Temma (IMakishima et al. Ill990bl ). Further more, its two higher 



harmonics have been discovered at energies of ^ 5 and ~ 75 keV (ICoburn et al. 1 12005 



Kreykenbohm et al. II2005I : iPottschmidt et al. II2005I ). Thus, X0331+53 is the second object 



which has three or more CRSFs in the spectrum, following 4U 0115+63 (jSantangelo et al. 



19991 : iHeindl et al. 1119991 ). However, compared with the results from 4U 0115+63, this object 
exhibited two intriguing differences (Mowlavi et al. 2006; Nakajima 2006; Tsygankov et al. 
2006, 2009). One is that the change in i?ai started at a higher luminosity, 1 x 10'^^ erg s^^. 



- 3 - 



at an assumed distance of 7 kpc fiNegueruela et al. Ill999l ). The other is that the second 



harmonic absorption feature persisted in the X-ray spectra through an outburst. 

Including the above cases of 4U 0115+63 and X0331+53, the behavior of the fundamen- 
tal CRSF has been extensively studied so far. In contrast, details of the second harmonic, 
including its possible luminosity dependence, are less understood, primarily because of lower 
data quality at higher energies. With this in mind, we analyzed the whole RXTE (ASM, 
PCA and HEXTE) data of X0331+53, acquired in a 2004-2005 outburst. As a result, we 
have discovered that the second CRSF energy Eia possibly depends more weakly on the 
luminosity than i?ai. The Ea2/Eai ratio was ~ 2.2 when the source was luminous, while it 
approached the nominal value of 2.0 toward lower luminosities. These results suggest that 
the two resonances take place at different heights in the same accretion columns, and that 
Es,2, formed at a lower height, provides a more reliable estimate of the surface field strength. 



OBSERVATIONS and DATA REDUCTION 



In order to study the luminosity-dependent changes of i?ai and i?a2 of X0331+53, here 
we utilized all the RXTE data, sets of this transient pulsar acquired in the 2004-2005 outburst. 
Figure [1^ shows the whole light curve of this outburst acquired with the All Sky Monitor 
(ASM; Levine et al. 1996) onboa rd RXTE. An abrupt brightening was detected on 2004 



Nover nber 25 ([Swank et al. II2004J ). and the 2 — 12 keV flux continued to increase up to 1.1 



Crab (IRemillard 1120041 ). As revealed by the ASM, the X-ray intensity declined toward middle 
of 2005 February, and exhibited a small recovery at the beginning of 2005 March, presumably 
synchronized with the 34.25 day orbital period ( IZhang et al. II2005I . and references therein). 
At the outburst peak, a number of target-of-opportunity observations were performed by 
INTEGRAL, Swift, and RXTE. 

From 2004 November 27 through 2005 March 27, 108 pointing observations were made 
with the Proportional Counter Array (PCA; Jahoda et al. 2006) and the High Energy X- 
ray Timing Experiment (HEXTE; Rothschild et al. 1998) on board RXTE. Based on the 
inspection of the operation status of the five proportional counter units (PCU0~PCU4), we 
have selected the data sets acquired with PCU2, which worked throughout this outburst. 
As reported by Pottschmidt et al. (2005), HEXTE cluster A did not function from 2004 
December 13 to 2005 January 14. Thus, we utilize only the data acquired with HEXTE 
cluster B. 



Using part of the RXTE data sets that were ac quired in the descent phase, the behav- 
ior of the fundamental CRSF was already studied ( iMowlavi et al. I l2006l : iNakajima I l2006l : 
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Tsygankov et al. II2006I ). In the present paper, we complement these studies of the descent 
phase, by analyzing those datasets which were left unused. In addition, we analyze the 
ascent-phase datasets to examine whether the CRS F energy exhibits any h ysteresis effects 
similar to those suggested by the 4U 0115+63 data (INakajima et al. Il2006a( ). 



During the 4 months of observation, several energetic background events, such as large 
solar flares and precipitations of high-energy particles trapped in the radiation belts, affected 
some of the data sets. In order to exclude such low-quality data, we inspected the electron- 
rate light curves, and have selected 86 data sets with good quality as listed in Table [U To 
analyze them, we utilized the HEAsoft version 6.0.5. We assume 1% systematic error for all 
energy bins of the PCA data (same as in Nakajima et al. 2006a). All of the errors presented 
in this paper are 90% confidence levels. 



DATA ANALYSIS AND RESULTS 



3.1. Analysis of Representative Spectra 



As shown in Figure [H the intensity of X0331+53 changed largely in the ASM, PCA 
and HEXTE energy bands during this outburst. In addition, the hardness-ratio exhibited 
an anti-correlation to the source count rates as shown in Figure [T]i. In order to investigate 
the luminosity related CRSFs changes through this outburst, we need to establish a uni- 
fied spectral modeling which can be applied consistently to the data sets at all luminosity 
levels. So, we first examined pulse-phase-averaged PCA and HEXTE spectra from several 
representative data sets. 

Since many of the data sets have rather short exposure, we selected four representative 
ones sets with relatively long exposures, denoted in Tabled] as Dec 2a, Dec 24, Jan 20, and 
Feb 13b; these represent the beginning, the peak, the declining phase, and the end of the 
outburst, respectively. Figure |2^ shows the background-subtracted PCU2 (3 — 20 keV) and 
HEXTE B (20 — 80 keV) spectra obtained on these four occas ions. The PCA back ground 



spectra were all estimated with the bright-background model (jjahoda et al. II2006I ). while 



the HEXTE backgroun ds were extracted from t he off-source position which is offset from the 
source by 1.5 degrees (IRothschild et al. 1119981 ). The four spectra, even in their raw forms, 
clearly reveal the prominent fundamental CRSF at ~ 30 keV, which was observed in the 
previous (Makishima et al. 1990a,b) and the present (Kreykenbohm et al. 2005; Pottschmidt 
et al. 2005; Mowlavi et al. 2006; Nakajima 2006; Tsygankov et al. 2006, 2009) outburst of this 
object. Furthermore, the brighter two spectra clearly show the second harmonic resonance 
at ~ 60 keV. 
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From early days (e.g. White et al. 1983), continuum spectra of accretion powered 
pulsars were approximated by a power-law modified by exponential cutoff. In this paper, we 



employ NPEX (Negatiy e and Positive power-law with Exponential) model ( iMihara Ill995 
Makishima et al. Ill999l ). which is described as 



N{E) = + A2E+^-') exp ( 



(1) 



where E is the X-ray energy in units of keV, Ai and ai are the normalization and photon 
index of the negative power-law, respectively, A2 is the normalization of the positive power- 
law, and kT represents the cutoff-energy in units of keV. This model successfully reproduced 
the continuum spectra of 4U 0115+63 (Mihara et al. 2004; Nakajima 2006; Nakajima et al. 
2006ab) at various luminosity levels. As already reported by several authors (e.g. Makishima 
et al. 1999; Coburn et al. 2002), the determinations of the cyclotron line parameters are 
affected by the continuum and the line modeli ng. Employing this mo del, we can compare 
the results of this work with the previous ones (INakajima et al. Il2006al ). Therefore, we have 
selected this model. 

We first attempted to fit the 4 spectra with the NPEX model, with Ai, A2, ai, and 
kT all left free. Since the PCA background estimations are only good to a few percent, 
we further adjusted the background normalizations down to a level of ~1%, following the 
method described in the RXTE cook hook^. However, as expected, the NPEX model gave 
acceptable fits to none of the four spectra. As revealed clearly by the data-to-model ratio 
in Figure [2]d, all the four spectra exhibit the fundamental CRSF strongly at ~ 30 keV. The 
second harmonic CRSF is observed not only in the brighter two spectra, but also in the two 
fainter ones. In addition, the spectrum on Dec 24 exhibits evidence of the third harmonic 
CRSF at ~ 75 keV, as already reported by Coburn et al. (2005), Kreykenbohm et al. (2005), 
and Pottschmidt et al. (2005). Below, we concentrate on the fundamental and the second 
harmonic. The fluo rescent Fe K^, lines at 6.4 keV with a width of ape = 0.5 k eV are also 



taken into account (jPottschmidt et al. II2005I : iNakajima II2006I : iTsygankov et al. II2006I ). 



In order to evaluate the fundamental and second CRSF parameters, we next intro- 



duce a cyclotron absorption moc 


el. Although some authors usee 


Gaussian-shaped a 


jsorp- 


tion cross section ( Co 


3urn et al. 


2002:lKrevkenbohm et al. 2004 


Pottschmidt et al. 


2005; 


Krevkenbohm et al. 


2005: Mowlavi et al. 2006: Kretschmar et al. 2006: Klochkov et al. 


2008) to reproduce the CRSF. we emplov the cvclotron absorption fCYAB) factor which 



has been used successfully (Clark et al. 1990; Makishima et al. 1990a, 1999; Mihara 1995; 
Mihara et al. 1998; Nakajima 2006; Nakajima et al. 2006ab; Tsygankov et al. 2006, 2007, 



^See [http://heasarc.gsfc.nasa.gov/docs/xte/recipes/cook_book.htinl 
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2009). The CYAB factor is described as 



Ci{E) = exp 



A {WE/Ea,Y 
{E - Eai)^ + 



1,2), 



(2) 



where Eai is the resonance energy, Wi is the resonance width, Di is the resonance depth, and 
i is the harmonic number (with i = 1 fundamental). Using high-quahty Suzaku data, Enoto 
et al. (2008) confirmed that eqJH can successfully reproduce the 36 keV CRSF of the X-ray 
pulsar Her X-1, whereas the Gaussian-shaped absorption cross-section is less successful. 

Since we already know that X0331+53 has the multiple CRSFs, we attempted to fit 
the data with an NPEX model multiplied by two CYAB factors (hereafter NPEXxCYAB2 
model). The values of Ai and Es2 were both left to vary independently, rather than con- 
strained as E^2= 2-Eai- As a result, the fits were much improved compared to the NPEX 
fit, as shown by the data-to- model ratios in Figure |2]:. Due to the complex shape of the 
fundamental CRSFs, however, the NPEXxCYAB2 model is not yet fully successful, with 
yf, > 1.2. Th i s fact was already reported from previous studies ([Makishima et al. I Il990al : 



Mihara II1995I : iKrevkenbohm et al. II2005I : iPottschmidt et al. II2005I ) 



In order to better describe the fundamental CRSF, Kreykenbohm et al. (2005) and 
Pottschmidt et al. (2005) introduced an additional Gaussian absorption (GABS) factor with 
a similar energy but a different width, thus constructing a nested double- Gaussian absorption 
profile. Therefore, let us try to explain away the fit residuals employing an additional GABS 
factor [^, which is described as 



G{E) = exp 



Tga exp 



1 fE-E, 



ga 



ga 



(3) 



where i?ga is the line-center energy, cXga is the fine width, and Tga is the optical depth at 
the line center. Although a more sophisticated CRSF modeling has been introduced by 
Schonherr et al. (2007), we employ the CYAB and GABS models in the present work, so 
that the results can be directly compared with previous studies of the luminosity dependence 
of CRSFs (Makishima et al. 1999; Mihara et al. 2004; Nakajima et al. 2006). 

We fitted the four spectra by the NPEXxCYAB2xGABS model, namely the same 
NPEXxCYAB2 model but further multiplied by eqj3l in which the three GABS parameters 



C"ga5 Tgaj 



were left free but £^ga was given an initial value close to E^^i. Due to low 



^The original GABS model in XSPEC vll.3.2p had a bug related to the energy binnings as reported by 
Kitaguchi et al. (2007). In this paper, we utilized an updated version of GABS model with the bug fixed, 
provided by Keith Arnaud. 
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statistics in the higher energies, we fixed i?a2 at 2i?ai in fitting the Dec 2a and Feb 13b data, 
while we left i?ai and i?a2 both free and independent for the brighter two data sets. Then, as 
illustrated in Figure lUi, the four spectra have been described successfully with this model. 
The derived best-fit parameters are summarized in table [2J To see the configuration of the 
nested cyclotron models for the fundamental resonance, the spectrum observed on Dec 24 
is presented in Figure [3] in its vFy form. Thus, the CYAB factor explains the gross shape 
of the fundamental CRSF, while the remaining narrow core is represented by the additional 
GABS factor. The value of £'ga agrees, within a few percent, with Eai {1 + {W^ /Eg^ Y\ 



or ~ 1.2-5'ai in the pres ent case, where the CYAB factor becomes deepest (IMihara I Il995 



Makishima et al. 1119991 ). For comparison, Kreykenbohm et al. (2005) and Pottschmidt et 
al. (2005) applied nested two GABS factors to the fundamental CRSF and obtained two 
Gaussian centroids which differs by 10 — 20%. 

Since the continuum and the CRSF factor couple strongly, the 90%-confidence errors 
of individual fit parameters, given in Table 2, might be significantly underestimated. To 
examine this concern, we present in Figure H] confidence contours between several pairs 
of the model parameters, obtained from the Dec 24 data. Thus, kT and the resonance 
energies are almost uncorrelated. Although the NPEX ai and the CYAB Eg^i exhibit some 
correlations, the single-parameter 90% error ranges are confirmed to adequately represent 
the two-dimensional confidence ranges, without being under-estimated. 

From these spectral analyses, we have established the unified model, which can re- 
produce the complex shape of the X0331+53 spectra regardless of the source luminosity. 
Hereafter, we use this model, namely, NPEX xCYAB2x GABS. 

Although our analysis was performed with the HEAsoft version 6.0.5, the PCA response 
generator has been updated in the latest HEAsoft version 6.8, and a new PCA CALDB 
has been released on 2009 December 2. To examine possible effects of these software and 
calibration updates, we analyzed the Dec 24 data with the HEAsoft version 6.8 and the latest 
PCA CALDB. As a consequence, the cyclotron resonance energies changed only ~1.0%, 
which is within statistical errors. Thus, we retain our results obtained with the HEAsoft 
version 6.0.5. 

Just for comparison, we attempted to replace the CYAB modeling of the two CRSFs 
with that employing three GABS factors, with nested two for the fundamental and the other 
for the second harmonic. The three parameters of the three GABS factors were all left free. 
As a result, the model gave, for example, an acceptable fit (x^ ~ 1.16) to the Dec 24 data. 
The obtained CRSF parameters are, E^a = 27.2+°;^, = S.lOtlH, = lAQ^oit, Eu = 
25.6^°i a,, = 3.23l°:S, = 0.381°;^^, E^ = 50.2^°^ = 7.08l°|^, and ='l.49l°Jt 
Here, the subscript of la and 16 specify the parameters of the nested two GABS model for the 
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fundamental resonance. However, we retain our original modeling using the CYAB model, 
because it gives generally better fits: for example, the difference in is 28 (for = 51) in 
the case of the Dec 24 spectrum. 



3.2. Analysis of the Date-sorted Spectra 



We apphed the model established in §3.1l to all the daily-averaged spectra, and studied 
the time evolution of the CRSF parameters. Since some data sets have insufficient statistics, 
the GABS model parameters did not converge. Following the obtained results in the previous 
subsection, we therefore fixed E^^ = 1.2Eai in all cases. The other procedures of the spectral 



analyses are the same as in §3.11 

As a result of this analysis, our model has given acceptable fits to all the data acquired 
in the major outburst spanning 2004 November through 2005 February. While the normal- 
ization parameters of the NPEX model varied in correlation with the luminosity, the NPEX 
fcT, which is thought to give a measure of th e electron temperature in the emission region 



( iMakishima et al. I Il999l: iMihara et al. I 20041). stayed rather constant like in the previous 



results on 4U 0115-1-63 flNakaiima et al. Il2006al : iTsvgankov et al. 1120071 ). Although the fun- 
damental CRSF parameters are well determined, the second CRSF parameters in some data 
sets became unconstrained due to insufficient data statistics in higher energies. In such cases, 
we fixed -E'a2 at 2_E'ai- 

In contrast to the November-February outburst data, the remaining three data sets 
were observed in March during the minor outburst recovery (Figured]). Since the source was 
rather faint on these occasions, the three spectra were reproduced with the NPEX model 
multiplied by a single CYAB factor describing the fundamental CRSF. Therefore, we do not 
discuss these data sets acquired in 2005 March. 

Figure shows the derived values of -Eai, together with 90%-confidence errors, as a 
function of the 3 — 80 keV source luminosity L^ ^q. As already reported by Mowlavi et al. 
(2006), Nakajima (2006), and Tsygankov et al. (2006, 2009), E^i thus changed from ~ 22 to 
~ 27 keV as -^3,80 varied between 5.0 x 10'^'^ and 3.5 x 10^^ erg s~^. When the data points in 
|5^ are fitted by a linear function of ^3,80 having a form of 

Eal=i?i?{l + W^l} , 



(4) 



its two parameters (with 90% confidence errors) were obtained as E^^' = 26.5 ± 0.1 keV 
and Li = (—23.0 ± 1.2) x 10^^ erg s~^, with xl = 1-3 {jy = 81). In addition, we found no 
hysteresis effects in the variation of Es,i between the ascent and descent phases. This inference 
is consistent with a recent report by Tsygankov et al. (2009), who used the same data set. 
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Although the values of E^i derived h ere are slightly discrepant with those measured from 



the same outburst by other authors (IKrevkenbohm et al 



2005: Pottschmidt et al. 2005 



Mowlavi et al. II2006I : iTsygankov et al. Il2006[ ). the difference is within ~ 10%, and can be 
attributed to the different modelings of the fundamental CRSF and the continuum. 

As shown in Figure [5]d, the second resonance energy i?a2, when determined indepen- 
dently, showed a considerably weaker dependence on L380. When these data points are 
fitted by another linear function of -^3,80 as 



E,2 = Ej^^{l + Ls,so/L2} 



(5) 



we obtain e!^ 



53.5 ± 0.4 keV and L2 

xt = 1-8 {jy 

differs significantly between £^ai and i?a2 (i-e 



-47.8 ± 5.6) X 10^* erg s~\ together with 
1 2.0-^^1'*! while the luminosit 
IL2I > \Li\ beyond their errors) 



61). Thus, Ej^^ is consistent with 2.0-^^1'*! while the luminosity dependence 



As a result of the above two correlations, the E.^2/ E.^i ratio, shown in FigureEt, increases 
from ~ 2.0 to ~ 2.2 as L^^q changes from 1.0 x 10'^'^ erg s~^ to 3.5 x 10^^ erg s~^. Since 
these ratios give = 2.8 {v = 62) when fitted with a constant value, they are inconsistent 
with being constant. In contrast, the fit became satisfactory with xt — 1-2 {'^ = 61), when 
we employ yet another linear function as 



E,2/E,i = K{l + Ls,8o/L2i} 



(6) 

The two parameters were obtained as = 2.0 ±0.1 and L21 = (34 ± 5) x 10'^^ erg s~^. 
Thus, the value of K is consistent with 2.0, while L21 remains finite. In other words, the 
Ea2/Eg,i ratios exhibit a statistically significant dependence on the 3-80 keV luminosity. In 
§3.3 and §4.2, we discuss possible artifacts on these results introduced by our choice of the 
CRSF modeling. 



3.3. Analysis of the Flux-sorted Spectra 

In order to investigate the behavior of the second CRSF under better statistics, we next 
carried out the flux-sorted analysis which was already performed successfully on 4U 0115+63 
(Nakajima 2006; Nakajima et al. 2006a,b). Specifically, we sorted the PCU2 and HEXTE 
cluster B data into 8 flux levels in reference to Figure O and co- added those data which fall in 
the same flux range. The flux was calculated every 64 sec, so that the flux sorting can catch 
up with short-term intrinsic variations, but not affected by photon-counting statistics. In 
addition, the whole data were divided into outburst ascent and descent phases (referring to 
Figure[6]), to examine the CRSF energy changes for possible hysteresis effects. The procedure 
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of the spectral analysis is the same as in §3.2[ but Eg^. was left free, because of the improved 
statistics. 

Figure [7^ and[H^ show the flux-sorted spectra of the ascent- and descent-phases, respec- 
tively. We again applied the NPEXxCYAB2xGABS model to these data, and found that 
the model gives acceptable fits as shown in Figure [7)d and [8]d. 

As summarized in Table [3l this analysis allowed us to accurately determine the parame- 
ters of the fundamental and second CRSFs. Figure [Hti' and b show the values of i?ai and i?a2, 
respectively, against L3 go- Thus, i?ai changed from ~ 22 to ~ 27 keV, while -Ea2 from ~ 49 
to ~ 54 keV, thus reconfirming the essential features of Figure [5] with higher confidence. In 
fact, like in the case of the date-sorted spectra (§3.2), the fit with eqJDto Figure gave 
E^^ = 21 A ± 0.6 keV, Li = (-16.5 ± 3.3) x lO^^ erg s~\ and xl = 0.5 {u = 12), while 
that with eqElto Figure resuhed in E^^ = 53.3 ± 0.5 keV, L2 = (-45.7 ± 7.0) x lO^^ erg 
s~^, and = 1-5 (i^ = 12). The errors associated with these fit parameters are larger than 
those obtained in the date-sorted analysis, because here we left Eg^ free while we previously 
fixed it at 1.2£'ai. Figure [Ht reconfirms that the resonance energy ratio i?a2/-£'ai changed 
from ~ 2.2 to the nominal value ~ 2.0 as the source luminosity decreased. The fit with eqE] 
indicates K = 2.0 ± 0.1, L21 = (21.2 ± 7.0) x 10=^^ erg s-\ and xl = 0.3 (z/ = 12). 

Through the date-sorted and flux-sorted analyses, we confirmed that the luminosity- 
dependent changes in ii^a2; E^i and the i?a2/-£'ai ratio are statistically significant. Then, 
does this conclusion remain unaffected even if considering systematic errors? An immediate 
concern, namely spurious effects due to coupling between the continuum and CRSF param- 
eters, can be ruled out because the error estimates are appropriate as indicated by Figure IH 
Another problem is that the values of E^i and -E'a2 do not pick up the deepest position of 
the CRSFs. Accordingly, we examined the behavior of the energy E^j = E^j |1 4- (Wi/ E^.j)'^} 



(i = 1,2) at which the absorption is expected to become deepest (IMihara Ill995l ). Figure ITU] 
shows luminosity dependence of E^i- Thus, the values of Eg^i also exhibit a clear luminosity 
dependence. By fitting these data points with eqJD and eqj5l we obtained the characteristic 
luminosities as Li = (-27.0 ± 11.8) x 10^^ erg s'^ and £2 = (-44.3 ± 8.0) x 10^® erg s"^ 
We hence reconfirm that the second resonance energy depends less steeply on the luminosity 
than the fundamental energy, although the difference between Li and L2 is less significant 
than that between Li and L2, and the average E^2/Ea.i ratio, 1.8-2.0 is smaller than 2.0 in 
contrast to the case of the i?a2/-£'ai ratios. 

The largest systematic effect that may possibly affect the £'a2/-E'ai ratio would be the 
choice of the continuum and CRSF models. Since we have already tried in §3.11 replacing 
the CRSF model with the GABS model, here we examine the choice of the continuum. 
In order not to be affected by the CRSF modeling, we fitted the lowest-energy (3-13 keV) 
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and highest-energy (60-80 keV) ends of the flux-sorted fdl and fd6 spectra, using a single 
NPEX continuum, and then normalized the overall spectrum to the best-fit model that 
was determined by each spectrum. The results, shown in Figure [TTh . reveal the absorption 
features in a relatively model-free manner, with the continuum approximately removed. 
There, we have adjusted the energy scale and the ratio scale between the two spectra, 
so that the fundamental CRSFs best overlap. Thus, the fdl ratio spectrum exhibits the 
2nd harmonic at a significantly higher energy than the fd6 ratio. This result visualize the 
luminosity-dependent change in the i?a2/-£'ai ratio. 

To examine the overall behavior of En2, E^i, and their ratios under a different continuum 
model, we employed a power- law cutoff (PLCUT) model (Coburn et al. 2002), namely a 
power-law multiplied by an exponential cutoff factor. Figure [TTb represents the same analysis 
as Figure [TTk conducted with this continuum model. Thus, the result again visualize that fdl 
has a higher Ea2/Eai ratio than fd6. We further fitted the fdl and fd6 spectra by multiplying 
these PLCUT continua with the same absorption factors as used so far, and obtained the 
results summarized in Table HI Thus, the fit goodness is comparable to those obtained 
with the NPEX continuum, and the values of E^i and Es,2 remain unchanged within errors, 
although the luminosity-dependent change in the i?a2/-E'ai ratio became less significant. 

As our further confirmation of the luminosity-dependent changes in the CRSF energies 
and their ratios, we fitted the fdl and fd6 spectra simultaneously, with the NPEXxCYAB2 
model. While the model parameters were generally allowed to differ between the two spectra, 
we imposed 4 stepwise constraints on the resonance energies. In the first step, we required 
Es.1 to be the same between fdl and fd6, and constrained ii^a2 to be twice i?ai- However, 
the fit was not acceptable, with x'^/'^ = 2365/111. Next, we allowed fdl and fd6 to have 
different values of E^^i, but retained the 1:2 harmonic constraint on both spectra. This gave 
an improved fit with x^/^ = 398/110, implying that the two spectra have different resonance 
energies. Third, we allowed the two spectra to have different values of Eg^i, but required them 
to have the same E^2/E^i ratio. The fit was then improved to x^/v = 321/110, together 
with Ea2/Es,i= 2.11 ± 0.08 which deviate from the the nominal value 2.0. Finally, we left 
all the four resonance energies free, and obtained a improved fit with x'^/^ = 282/108. This 
means that the i?a2/-E'ai ratio is different not only from the harmonic condition, but also 
between the two spectra. (The fit is still unacceptable because we did not include the GABS 
factor to model the narrow core of the fundamental resonance.) 

Finally, we repeated the same 4-step analyses, but using the PLCUT continuum. Then, 
the fit goodness from the 4 steps were x'^/'^ = 2091/111, x^/^^ = 1119/110, x^/i^ = 341/110, 
and x^/ ^ = 254/108. The implication is essentially the same as that derived with the NPEX 
continuum. 
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3.4. Additional Analyses 

Since the flux-sorted spectra have rather good statistics, we examined them, up to 110 
keV, for evidence of the third harmonic CRSF which was already reported by Kreykenbohm 
et al. (2005) and Pottschmidt et al. (2005). For this purpose, we further multiplied the 
NPEXxCYAB2xGABS model by an additional CYAB factor centered at ~ 75 keV, and 
applied it to the flux sorted spectra. Then, two data sets, fal and fd2, gave statistical 
signiflcant (at ~90%) evidence for the third CRSF, while the fundamental and second- 
harmonic CRSF parameters remained essentially the same as before. As a representative 
case, the third resonance parameters derived from the fal data are -E'a3= 73^1 keV, = 6^1 
keV, and = 3^^^. Th ese parameters are comparable to the previous resu lts derived 



from the same RXTE data ( iKreykenbohm et al. II2005I : iPottschmidt et al. II2005I ). The third 



resonance energy was found to change no less than by ~ ±3%, as the 3 — 80 keV luminosity 
varied from 3.16 to 3.59 x 10^^ erg s~^. Over this luminosity range, the E^/ E^2 ratio stayed 
at an average value of 1.46 within a typical uncertainty of ±0.05. In contrast, the E^/E^i 
ratio is systematically higher, 3.2. Therefore, E^^ and i?a2 are consistent with being in the 
2:3 harmonic ratio, while i?ai is inferred to be systematically lower, and depend signiflcantly 
on the luminosity. 

In addition to the investigation of the resonance energy changes, we examined whether 
the depths of the cyclotron resonance, (i = 1, 2) and Tga, depend on the luminosity. Figure 
[T2] shows these depth parameters, against the 3 — 80 keV luminosity. Both the fundamental 
and second CRSF depths, Di and D2 respectively were found to decrease toward higher 
luminosities. Although Tga has large errors, it also exhibits the same tendency. These 
luminosity dependence in Di and D2 are statistically signiflcant, because fltting the Di 
and D2 measurements with constant values gave reduced chi-square of 2.11 {v = 13), 5.88 
= 13), respectively. 



4. DISCUSSION 

4.1. Modeling of the Spectra 

We analyzed the whole set of RXTE PCA and HEXTE data of X0331+53, covering 
the 2004—2005 outburst, with two objectives in mind. One is to examine whether the 
luminosity-anticorrelated changes in the fundamental CRSF energy, observed during the 
outburst descent phase (Mowlavi et al. 2006; Nakajima 2006; Tsygankov et al. 2006), is 
also present during the outburst rising phase. This agrees with the conclusion of Tsygankov 
et al. (2009). The other is to examine the reported second harmonic resonance of this 
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pulsar (Coburn et al. 2005; Kreykenbohm et al. 2005; Pottschmidt et al. 2005) for possible 
luminosity-dependent changes, and compare the results with those of the fundamental CRSF. 

Following our experience (e.g., Mihara 1995; Makishima et al. 1999), we adopted the 
NPEX model of eq{T]to represent the underlying 3—80 keV PCA and HEXTE continua, while 
two CYAB factors, eqJH to account for the fundamental and second harmonic cyclotron 
features. Furthermore, to express the com plex shape of the deep fundamental cyclotron 



featu r e as revealed in the prey ious studies (IMakishima et al. Ill990al : Kreykenbohm et al. 



20051 : iPottschmidt et al. Il2005[ l. we introduced a fine-tuning GABS factor, eqj3l which is 
nested with the fundamental CYAB factor. This NPEXxCYAB2xGABS model has given 
acceptable fits to all the daily-averaged and fiux-sorted spectra. Furthermore, some of the 
spectra exhibited statistically significant evidence for the third harmonic resonance, as al- 
ready reported by Coburn et al. (2005), Kreykenbohm et al. (2005), and Pottschmidt et al. 
(2005). 

As exemplified in Figure [12] and Table [3l our fits generally imply that the fundamental 
and second harmonic resonances both have optical depths of 1.0— 2.0, with comparable CYAB 
widths of ~ 8 keV. The fundamental CRSF is explained mainly by the CYAB model, with 
the additional GABS factor carrying ~ 30% of the optical depth when the source is luminous. 
Replacing the two CYAB models with two GABS models made the fits significantly worse 
(e.g. xl= 3.6 in fa2 spectrum fitting), mainly because the observed flux falls considerably 
below the NPEX continuum, even at the 35 — 45 keV range which is in between the two 
cyclotron troughs (see Figure [3]). This effect cannot be reproduced adequately by the GABS 
model which is symmetric between the red and blue sides, while it can be successfully 
accounted for by the CYAB factor of which the wing is more extended in the blue side due 
to the E"^ factor in eqlH 



4.2. Changes of the Cyclotron Resonance Energies 



By applying the NPEX xCYAB2x GABS model to the daily-averaged and fiux-sorted 
PCA+HEXTE spectra, we reconfirmed the reports by Mowlavi et al. (2006), Nakajima 
(2006), and Tsygankov et al. (2006, 2009), that the fundamental CRSF energy E^i decreases 
as the source gets more luminous (Figure|9]). Thus, X0331+53 becomes a second binary X-ray 
pulsar, after 4U 0115+63, of which i?ai correlates negatively with the source luminosity. Like 
the case of 4U 0115+63, this effect observed from X0331+53 can be explained by presuming 
that the cyclotron resonance "photosphere" gets higher in the accretion column, a.s the 



source luminosity inc r eases and hence the column becomes taller (IBasko fc Sunyaev. Ill976 



Burnard et al. Ill99ll : iMihara et al. Il2004| ). These cases make a contrast to the behavior 
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of Her X-1, in which E^a is reported to have been varyin g in a positive correlation with 



the luminosity (iGruber et al. I I2OOII : IStaubert et al. 1 120071 ) , at least while the source is in 
a sub-Eddington reginie. On the other hand, the recent works of Her X-1 (IStaubert et al. 



20071 : iKlochkov et al. 1 12008| ) reported that the anti-correlation behavior might exist when 



the source was in the super-Eddington regime. 

In 4U 0115+63, the luminosity dependence of £^ai showed some hysteresis effects between 
the ascent and descent phases of the outburst (Nakajima et al. 2006). In contrast, we found, 
in agreement with Tsygankov et al. (2009), no such effects from X0331+53; in this source, 
the value of i?ai can be considered as a single- valued function of the source luminosity. 

As the most important result of the present work, our NPEXxCYAB2xGABS modeling 
leads to an inference that the second resonance energy E^2 also correlates negatively with 
the source luminosity, but significantly more weakly than E'ai does. The two resonance 
energies decrease by A£'ai/£'ai ~ 12% and Ai?a2/-E'a2 ~ 5%, as the 3 — 80 keV luminosity 
increases from 1.0 x 10^^ to 3.5 x 10^^ erg s~^. Consequently, the £'a2/-E'ai ratio is consistent 
with the nominal value of 2.0 when the source is dim, while it increases to ~ 2.2 toward 
the outburst peak (Figure (Ht). These effects are somewhat model dependent, and become 
less significant when the absorption profile and/or continuum shape are modeled in different 
ways. Although these effects were already noticed by Tsygankov et al. (2006), their spectra 
had rather low statistics in higher energies, and hence the second CRSF parameters were 
not well constrained. 



While the relativistic effects (lAraya fc Harding. Ill999[ ) predicts the i?a2/£'ai ratios to be 
lower than the harmonic value of 2, actual measurements do not necessarily agree. The ratio 
exceeds 2.0 in some cases, and fall in the range 2.1 ~ 2.2; for example, 4U 1907+09 showed 
Eai ~ 18 keV and fi^ ^ 38 keV dCusumano et al. Ill998h: Vela X-1 s howed E^i ~ 24 keV 
and Ea2 ~ 52 keV jMakishima et al. Ill999l: IXrevkenbohm et al. II2OO2I ): and A0535+26 E^i 



50 keV and i?a2 ~ 110 keV ( Kretschmar et al. 20061 ). In contrast, the ratios lower than 



2.0 are reported in other cases, including the present RXTE data of X0331+53 themselves 
analyzed by some other authors (e.g., Pottschmidt et al. 2006) who generally took the deepest 
positions of the CRSFs as the resonance energies 



A dominant origin of the reported scatter in the E.^2 / -E'ai ratio is presumably the use of 
different (e.g., CYAB vs. GABS) modeling of the CRSF profile. In this sense, our result of 
Ea2/Eai> 2 is not free from this systematic problem. In fact, we find -E'a2/-E'ai<2 if instead 
using Es2/Es_i. Nevertheless, our result has some pieces of supporting evidence. One is 
that the CYAB modeling is more successful than that with the GABS factor ( §3.ip . when 
combined with the NPEX continuum. Another is that the third to second resonance energy 
ratio, Es,3/En2, is consistent with 1.5, while E^^s/Es,! is higher than the nominal value of 3.0. 
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Yet another support is provided by a physically self-consistent interpretation to be described 
in the next subsection. Even admitting that the absolute values of i?ai and might be 
subject to the systematic modeling uncertainly, the luminosity- dependent changes in the 
resonance energy ratio is considered to be robust as already shown in §3.31 



4.3. A Possible Explanation of the Observed Effect 

Let us consider a possible explanation of the luminosity-dependent changes in E^i, -E'a2, 
and their ratios. According to Mihara et al. (2004) and Nakajima et al. (2006a), the observed 
changes in E^^i of 4U 0115+63 can be explained in terms of those in the height of cyclotron 
resonance "photosphere" , in combination with a dipole field geometry. In this scenario, the 
photosphere height above the neutron star surface, hr, can be estimated as 



K fE, 



-Rns V Est 



-1/3 



- 1 , (7) 



where i?NS is the neutron star radius, Egi is the resonance energy to be observed on the 
neutron-star surface, and i is the harmonic number. 

Substituting the E'ai and £'a2 values derived from the fiux-sorted analysis ( §3.3p into 
eqJTJ hr of the fundamental and second resonances have been estimated as shown in Figure 
[131 Here, we employed E^i = 27 keV and Es2 = 54 keV, referring to the lowest- luminosity 
state in Figure |9l Also a typical value of i?NS = 10 km was employed. By translating 
the measurements into the photospheric heights, the figure thus yields two important im- 
plications. One is that the photosphere of the fundamental CRSF gets higher as the source 
luminosity increases, up to 800 m, or ~ 8% of -Rns- The other is that the second resonance 
photosphere, though increasing with the luminosity as well, is located closer to the surface 
than that of the fundamental, just reflecting the deviation of the Eg,2/E^i ratio from 2.0. 
The height difference between the two photospheres reach ~500 m when the source is most 
luminous. Like the estimation of hr described above, we also calculated the resonance heights 
using the asymptotic energies which are derived from E^i. Assume E^^i = 31.0 keV and 
Es2 = 62.0 keV, the luminosity- dependent height variations turned out to be ~ 500 m for 
Eai and ~ 400 m for E^2- 

The altitude difference between the two photospheres can be explained qualitatively in 



the following manner. Theoretically ( iHarding fc Daugherty. Ill99ll : lAraya &: Harding. Ill999 



and references therein), the fundamental CRSF is predicted to have a cross section which 
is ~ 10 times larger than that of the second resonance. Then, we expect the fundamental 
photosphere to be formed nearly at the top of the accretion column. In contrast, the second 
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harmonic photosphere will appear at lower altitudes, where the column density integrated 
along our line of sight is expected to become higher. These effects have already predicted 
theoretically by Nishimura (2008). Then, a detailed modeling of the ii^a2/-£'ai ratio may 
provide a valuable probe into the density distribution along the accretion column, and its 
dependence on the mass accretion rate, although such a work is beyond the scope of the 
present paper. 

Assuming that the luminosity- dependent changes in the resonance energies are caused 
by the photospheric height variations, and not, e.g., by relativistic effects ( §4.2p . the above 
physical picture further gives an a posteriori support to our finding of Eia/Eg^i> 2.0. In 
fact, if this ratio fell below 2.0 toward lower luminosities, we would have to conclude that 
the 2nd resonance photosphere is located higher than that of the fundamental, leading to a 
physically unrealistic condition. 



4.4. Comparison with 4U 0115+63 



Following our prey ious results on 4U 0115+63 (IMihara et al. I |2004J : iNakajima 1 12006 



Nakajima et al. Il2006al ). the present work provides valuable information as to the luminos- 



ity dependent changes of the CRSF parameters in accreting X-ray pulsars. Although the 
fundamental CRSF of X0331+53 required the additional GABS factor (which was not needed 
in 4U 0115+63), the CYAB and GABS centroids varied similarly as shown in Figure [9] and 
[T2l We therefore consider that this slight difference in the CRSF modeling does not hamper 
a direct comparison of the results from the two sources. 

Let us fist compare the luminosity dependence of E^i in the two objects. As shown in 
Figure [H^, the resonance energy of X0331+53 depends on the luminosity much less steeply 
than that of 4U 0115+63, even though the sense of dependence is the same between them. 
Specifically, the dependence is expressed as {AEai/ Eai) / {ALx / Lx) ~ 0.08 in X0331+53, 
while {AEai/Eai)/iALx/Lx) ~ 0.26 in 4U 0115+63. 

Assuming that the change in i?ai is caused by variations in the height of cyclotron 
resonance "photosphere", the above difference suggests that the accretion column height 
of 4U 0115+63 responds more sensitively to changes in the mass accretion rate. This 
inference is reinforced by Figure UM ). which directly relates the estimated photosphere 
heights with the luminosity. This difference, in turn, may be attributed to differences in 
the accretion column shape. Indeed, as shown in Figure [151 the resonance depths D of 
X0331+53 depend negatively on the resonance w idth W/Ea,, while the correlation is opposite 



in 4U 0115+63. According to previous studies ( llsenberg et al. Ill998l ; iKreykenbohm et al. 
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2004 : iNakajima et al. Il2006ai and references therein), this is considered to indicate that the 
accretion column in 4U 0115+63 has a tall cylindrical shape, whereas that in X0331+53 a 
fiat coin-like shape. 

The two sources are inferred to differ not only in the shape of the accretion col- 
umn, but also in its area Ac,r,\. We e xpect the observed luminosity to be expressed as 
Lx = AcoiCsbT^q (jBildsten et al. 1119971 ). where ctsb is the Stefan-Boltzmann constant, and 



Tefj i s the accretion column t e mperature which can be approximated by the observed NPEX 



kT (IMakishima et al. I Il999l : Mihara et al. 1 120041 : iNakajima et al. Il2006ah . By substitut- 



ing the derived parameters into this relation, we actually find X0331+53 to have a 1.43 
times larger Am] than 4U 0115+63. The relatively shallow spin modulation in X0331+53 
( iPottschmidt et al. Il2005[ ) supports this idea. Therefore, the accretion column in X0331+53 
is considered to have a coin-like shape and a rather large area, so that its height depends 
only weakly on the mass accretion rate. 

Let us finally discuss the behavior of the resonance depths. As already mentioned in the 
end of section IX^ Di, D2, and Tga of X0331+53 have all been found t o decrease toward higher 
luminosity. Thi s change is in the sam e sense as those in 4U 0115+63 (INakajima et al. Il2006al ) 
and GX301-2 ( lOkada et al. 1 120041 ). One possible mechanism causing such a dependence 
may be related to the appearance of higher harm onic resonances. As argued by several 



authors ( INakajima 1 12006 



Nakajima et al. I l2006al ) , an increased luminosity may give rise 



to significant higher harmonic resonances, which t end to make the fundaraental resonance 
shallower through so-called "two-photon" effects ( lAlexander fc Meszaros. Ill99l[ ). In fact. 



the third CRSF of X0331+53 was confirmed to appear when the source becomes luminous. 

In summary, we have quantified luminosity dependence of the fundamental and second 
resonance energies, and discovered clear luminosity-dependent changes in their ratio. In 
addition, we have discovered that the behavior of X0331+53 is qualitatively similar to that 
of 4U 0115+53, but differs quantitatively. 
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Table 1. The log of RXTE observations of X0331+53 in the 2004-2005 outburst. 







PCA 


HEXTE 








PCA 


HEXTE 


Date 


Start/End Time* 


Rate*" 


Rate'= 


Date 


Start /End Time* 


Rate** 


Rate" 


(2004/2005) 


(UT) 


[c s-1] 


[c s-1] 


(2004/2005) 


(UT) 


[c s-1] 


[c s-1] 


Nov 27 


14:06/14:13 


541 ± 2 


35.9 ± 1.4 


Jan 


lid 


14:48/15:52 


2049 ± 4 


73.3 ± 0.4 


Dec 1 


06:54/07:41 


848 ± 2 


44.5 ± 0.4 


Jan 


14 


08:43/15:59 


1841 ± 4 


66.8 ± 0.2 


Dec 2a 


03:20/04:12 


973 ± 2 


49.9 ± 0.4 


Jan 


15a 


00:28/04:33 


1789 ± 4 


64.9 ± 0.2 


Dec 2b 


19:16/20:11 


1042 ± 2 


52.2 ± 0.4 


Jan 


15b 


08:19/15:14 


1770 ± 3 


64.4 ± 0.2 


Dec 3 


18:52/00:25 


1121 ± 2 


54.1 ± 0.2 


Jan 


15c 


20:57/21:17 


1708 ± 4 


61.7 ± 0.5 


Dec 4a 


01:00/01:30 


1100 ± 2 


52.1 ± 0.5 


Jan 


15d 


22:28/03:59 


1707 ± 3 


61.3 ± 0.2 


Dec 4b 


02:34/03:14 


1141 ± 2 


55.2 ± 0.4 


Jan 


16a 


09:30/14:29 


1623 ± 3 


60.7 ± 0.2 


Dec 4c 


05:52/06:25 


1155 ± 2 


54.3 ± 0.6 


Jan 


16b 


20:33/02:12 


1604 ± 3 


59.5 ± 0.2 


Dec 4d 


09:15/09:39 


1189 ± 2 


57.2 ± 0.5 


Jan 


17a 


02:49/03:36 


1578 ± 3 


59.2 ± 0.4 


Dec 4e 


12:00/18:02 


1176 ± 2 


56.6 ± 0.2 


Jan 


17b 


07:31/13:00 


1541 ± 3 


58.1 ± 0.2 


Dec 4f 


18:21/00:00 


1204 ± 2 


56.3 ± 0.2 


Jan 


17c 


20:10/00:13 


1608 ± 3 


56.8 ± 0.2 


Dec 5a 


00:36/01:08 


1219 ± 2 


56.4 ± 0.5 


Jan 


18a 


00:50/01:55 


1475 ± 3 


56.1 ± 0.3 


Doc 5b 


13:11/14:16 


1287 ± 3 


60.0 ± 0.4 


Jan 


18b 


07:08/12:58 


1316 ± 3 


56.0 ± 0.2 


Dec 5c 


17:55/19:15 


1328 ± 3 


60.7 ± 0.5 


Jan 


18c 


21:17/01:24 


1433 ± 3 


54.8 ± 0.2 


Dec 6 


17:36/18:51 


1443 ± 3 


65.7 ± 0.5 


Jan 


19a 


02:01/02:46 


1440 ± 3 


55.0 ± 0.4 




1 / lUO/ 1 / lOl 




oy. ( i u.o 


Jan 


19b 


n^- /nsi-A 1 

\JO.Oi) / UoKl-l. 






Dec 13 


18:30/19:11 


2392 ± 5 


88.9 ± 0.6 


Jan 


20 


07:53/11:30 


1366 ± 3 


54.9 ± 0.2 


Dec 14 


17:25/18:21 


2557 ± 5 


93.7 ± 0.4 


Jan 


21 


07:42/08:21 


1295 ± 3 


53.1 ± 0.4 


Dec 15a 


12:18/14:20 


2609 ± 5 


95.7 ± 0.3 


Jan 


23 


05:07/06:32 


1201 ± 2 


50.3 ± 0.3 


Dec 15b 


14:51/19:42 


2614 ± 5 


96.9 ± 0.2 


Jan 


24a 


06:23/13:49 


1183 ± 2 


49.5 ± 0.2 


Dec 15c 


20:36/21:17 


2604 ± 5 


95.4 ± 0.5 


Jan 


24b 


14:13/15:07 


1165 ± 2 


48.9 ± 0.4 


Dec 16 


02:26/03:02 


2665 ± 5 


96.2 ± 0.5 


Jan 


25 


07:26/10:36 


1099 ± 2 


48.4 ± 0.2 


Dec 17 


11:44/12:51 


2789 ± 6 


97.2 ± 0.4 


Jan 


28 


10:57/11:43 


943 ± 2 


43.7 ±0.4 


Dec 18 


01:41/02:53 


2851 ± 6 


100.5 ± 0.5 


Jan 


29 


14:25/14:51 


927 ± 2 


43.9 ± 0.6 


Dec 19a 


14:14/15:10 


2870 ± 6 


100.9 ± 0.4 


Jan 


31 


08:48/09:13 


858 ± 2 


41.8 ± 0.5 


Dec 19b 


20:08/20:32 


2859 ± 6 


98.1 ± 0.7 


Feb 


2 


08:09/08:23 


780 ± 2 


37.8 ± 0.7 


Doc 20 


23:58/00:58 


3017 ± 6 


105.4 ± 0.5 


Feb 


4 


05:39/05:54 


780 ± 2 


38.5 ± 0.6 


Doc 22 


02:49/03:32 


3023 ± 6 


103.9 ± 0.7 


Fob 


6 a 


13:34/16:19 


681 ± 1 


35.2 ± 0.2 


Dec 24 


09:20/16:11 


3163 ± 6 


101.7 ± 0.2 


Feb 


6b 


18:16/18:45 


696 ± 1 


35.8 ± 0.4 


Dec 26 


08:32/15:48 


3120 ± 6 


100.0 ± 0.1 


Feb 


6c 


23:00/23:36 


665 ± 1 


35.7 ± 0.4 


Dec 28 


05:28/06:21 


3119 ± 6 


99.0 ± 0.5 


Feb 


8 


06:46/07:34 


618 ± 1 


32.9 ± 0.3 


Dec 29a 


05:52/06:30 


3078 ± 6 


99.5 ± 0.4 


Feb 


9 


18:02/18:14 


545 ± 2 


30.8 ± 1.0 


Doc 29b 


21:52/22:29 


3098 ± 6 


102.4 ± 0.6 


Fob 


10 


14:52/15:26 


516 ± 1 


29.6 ± 0.5 


Doc 3D 


13:35/14:20 


3040 ± 6 


97.7 ± 0.5 


Feb 


12a 


01:23/04:10 


478 ± 1 


24.2 ± 1.1 


Jan (ia 


06:23/06:38 


2582 ± 5 


85.5 ±0.8 


Feb 


12b 


17:04/21:06 


446 ± 1 


24.4 ± 0.3 


Jan 6b 


15:15/16:04 


2626 ± 5 


87.4 ± 0.4 


Feb 


13a 


00:58/01:45 


440 ± 1 


25.6 ± 0.6 


Jan 8a 


10:11/10:34 


2475 ± 5 


79.0 ± 0.6 


Feb 


13b 


04:12/08:18 


427 ± 1 


25.9 ± 0.2 


Jan 8b 


14:25/14:56 


2490 ± 5 


83.7 ± 0.6 


Feb 


13c 


11:45/14:10 


401 ± 1 


24.4 ± 0.2 


Jan 8c 


17:41/18:25 


2445 ± 5 


83.3 ± 0.5 


Feb 


13d 


18:01/19:01 


390 ± 1 


23.6 ± 0.3 


Jan 10 


10:57/16:02 


2157 ± 4 


76.5 ± 0.2 


Feb 


15 


12:07/14:56 


329 ± 1 


21.0 ± 0.3 


Jan 11a 


01:01/01:24 


2112 ± 4 


72.7 ± 0.5 


Mar 7 


07:07/08:06 


33.8 ± 0.1 


2.6 ± 0.2 


Jan lib 


08:20/09:18 


2159 ±4 


74.2 ± 0.4 


Mar 9 


07:54/08:10 


75.9 ± 0.5 


5.1 ± 0.7 


Jan 11c 


11:32/12:32 


2106 ± 4 


74.0 ± 0.4 


Mar 


18 


01:34/01:50 


12.4 ± 0.2 


0.8 ± 0.4 



^Start and end time (UT) of the PCA observations. 

^Background subtracted count rates of PCU2 in the 3-20 keV energy range. 
•^Background subtracted count rates of HEXTE cluster B in the 20-80 keV energy range. 
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Table 2. The best-fit parameters of NPEXxCYAB2xGABS model, determined by the 
four representative spectra shown in Figure [2l 



Obs. Date 



parameters Dec 2a Dec 24 Jan 20 Feb 13b 



1" 

iron 


1.54 ±0.26 


"•"^-0.99 


2 c;q+o.36 


65+°-ii 


ai 


-0.22 ±0.07 


"-"■^^-o.io 


-0 so+°-^° 


07+0-09 


kT (keV) 


6.25+?:«^ 


5 62+°-^'^ 


c- C-7+O.9O 
-0.79 


6 16+°-^^ 
"•-■-"-o.ei 


Eai (keV) 


25.8+?i 


21.1+^-? 


24.4+°| 




Wi (keV) 




10.71^?j? 


8.611?-°° 


9.83l}i^ 


Di 


1 00+0.27 

-'-■'-'^-0.29 






1 05+°-^^ 

i.UO_o.27 


Ea2 (keV) 


51.7 (fixed) 


49.5 ±0.3 




50.8 (fixed) 


W2 (keV) 


7.78i6;i9 


4.86^}:^? 




9.83(fixed) 


D2 


1 77+0.75 

-0.90 


1.35 ±0.12 


1 qc:+6.99 


1 qn+0.47 

-'-•^"-'-0.30 


^ga (keV) 


29.8+J-° 


26.4 ±0.3 


90 Q+0-5 

Z, J.O_g 4 


90 7+0-5 
' -0.3 


o-ga (keV) 


2.4il?:i 




2 70+°-^2 


9 /19+0.56 
-^•^^-0.43 


Tga 


o.9ot°:i 


0.48l°i^ 




1 qq+0.46 




1.37 


3.65 


1.76 


0.65 




0.77(52) 


0.55(51) 


0.66(51) 


0.67(52) 



'^In units of 10 ^ photons cm ^ s ^. 

''In units of xlO^^ ergs s~^ in 3 — 80 keV. 



Table 3. The best-fit NPEXxCYAB2xGABS model parameters, determined by the flux-sorted spectra. 



NPEX CYAB GABS 

flux ST kT ~i Wl dI W2 7^ 7^ Lx° 

level (keV) (keV) (keV) (keV) (keV) (keV) (keV) 



ascent phase 




fal 


-0.23+°- 


13 

.11 


fa2 


-0.24i°: 


12 
.10 


fa3 


-0.26+°- 


12 
.09 


faS 


-0.24+°; 


13 
.10 


fa6 


-0.20+°' 


10 
.09 


fa7 


-o-2si°: 


11 
.06 



6.43 
5.46 
6.03 
6.58 
6.60 
6.52 



+0.14 
0.12 
+0.17 
0.15 
+0.50 
0.39 
+0.71 
0.50 
+0.25 
0.23 
+5.14 
1.36 



20, 


.9+_l 


1 

.0 


10.791111 







13 
.16 


49.3 ± 0.3 


4, 


.33lS 


88 
.86 


, ,, + 0,10 
^■^^-0,09 


26. 


4+0 


4 
.3 


3, 


,39lS; 


46 
,48 


55+°-2^ 
u,oo_Q 21 


3, 


,69 


0, 


,81 


22. 




1 
.9 




1 


.021°: 


20 
.21 


49.5 ± 0.3 


4, 


.40+1 


18 
.10 


1 43+°'l= 
^^•'-0.12 


27. 


4+0, 
*-0 


.5 
.3 


3, 




49 
,56 


o.68lg;^,* 


3, 


,26 


0, 


,69 


23. 


.otl 


.0 
.2 


io.2ittll 


1 


.221°; 


22 

.29 


50.3l°;J 


6, 


.05+1 


59 

.08 


1 qq+029 
-^■^^ — 0.27 


27. 


6+° 
■°-o 


.5 
.3 


3, 




78 

,83 


0-67lg;re 


3, 


,07 


0, 


,71 


23, 


r + 1 
6_2 


2 
.0 


a 00 + 2. 64 
a.sj_2.oo 





.9612: 


24 
31 




5, 


,35+t 


65 
,27 


-^■°'^-0,41 


29. 


1+0, 


.4 
.3 


3, 


,23lS; 


59 
,58 


1 03+°'"' 


1, 


,79 


0, 


,67 


24. 


.3+° 


7 

.1 


^■■^'-1,29 


1 


■04lS: 


18 
16 


61.1 ± 0.4 


4, 


,63+i 


47 
,29 


1 aa+°-3i 

l'****-0,23 


29. 


„+0, 
■■5_o 


.3 
.2 


2, 


.93lS: 


34 
,41 


0.98 ± 0.27 


1, 


,63 


0, 


,66 


26. 




.0 
.8 




1 


.621°; 


33 

.44 


5i.8l|;^ 


5, 


.76l^' 


,22 
,75 


3.63l«°,™ 


31. 


■all 


.1 
.3 


0, 


,9ilg; 


60 
,42 


9.16(fixed) 


1 


,10 


0, 


,64 



descent phase 



fdl 


n 00+0.13 


6, 


,561°; 


51 

,36 


21, 


etl 


2 
,3 


I0.00l^;g| 


1 


.031° 


20 
.37 


49.3 ± 0.9 


6, 


nil 


70 
,95 


1 09+°'^° 


26.51°;^, 


q qO + 0.95 

■^.■J^ — 0.95 


o.5ol°:|g 


3, 


,54 


0, 


,64 


fd2 


„ ,Q+0.12 
"••^''_0.09 


6, 




84 
,72 


22, 


.9+_l 


,0 
,7 


8.30+f™ 


1 


.121° 


37 
.56 


49.7+i« 


8, 


,61+? 


08 
,19 


, „„ + 1.23 
l-^'^-0.57 


27.8ll;i 


3.38ll;°/7 


o.6ilg:fo 


3, 


,16 


0, 


,63 


fd3 


-0 24+°'^° 
"■■^*-0.09 


6, 


.09lg; 


32 
,28 


22, 


.9+1 


,8 
,1 


8.37+l-« 


0, 


.991° 


.21 
.25 


Bo.3i°;i 


3, 


40+1 


,39 
,40 


1 61+2™ 
^•"'^-0.34 


28.3l°;J 


2.92l°;« 


n 01+0.33 

0.01-0.28 


2, 


,72 


0, 


,80 


fd4 


-0 26+°'l° 
"•■^0-0.08 


6, 


,37lg; 


32 
,27 


23, 


q+0 


,7 
,9 


9.03ll;f, 


1 


.181° 


.18 
.17 


61.7 ± 0.5 


4, 


,94+1 


,79 
,68 


, ofi+0.23 
l.»O_0.23 


28.61°;^ 


„ __+0.3S 
'-0.51 


0.73 ± 0.26 


2, 


,23 


0, 


,61 


fd6 


-0 26+° " 
"••^0-0.08 


6, 


,3ilg; 


23 
,25 


23, 


.71?' 


,6 
,0 


^■^■'-1.05 


1 


.011° 


14 
.16 


51 6+°'* 
°^-°-0.5 


4, 


,56+1 


33 
,52 


, 00 + 0. 24 
i.»»_0 22 


28.9 ± 0.2 


2.94 ± 0.28 


„ qt.+0.23 
".^°-0.19 


1, 


.87 


0, 


,68 


fd6 


-0 3l+° 07 
"■■^^-0.06 


6, 


,28lJ; 


04 

,55 


25, 


9 + 0, 
"'-0 


6 
,8 


8 99+'^ " 
"■^"-i.oi 


1 


.35l°„ 


19 
,18 


52.41J;? 


7, 


.95+2 


62 
,29 


90 + 0,56 


29.3 ± 0.3 


„ 97+0.55 
^■'"-0.74 


n RQ+0.3° 

0-89_0.40 


1, 


.46 


0, 


,44 


fd7 


-0 04+° °'' 
"■^^-0.10 


6, 


■54lS: 


70 
,37 


25, 


7+0, 


5 
,9 


'■f«-0,97 


1 


.16lS 


22 
,26 


51.81J;? 


5, 


.74+_l 


62 
,32 


1 96+°" 


30.il°;J 


^.OD-0.58 


97+°.'lE> 
^.^'-0.34 


0, 


.79 


1, 


,08 


fd8 


07+° °'' 
-0.06 


7, 


■24l?: 


88 
,22 


26, 


0+0, 

■»-o 


,5 
,9 


10.031?;°^ 


1 


■^r+_° 


,36 
,30 


54.4l-.i 


8, 


,241° 


,40 
,15 


4 09+'*'** 
^■"''-1.34 


29.2 ± 0.7 


l-'°-0.97 


1 24+1-^1 
^■''*-0.93 


0, 


.66 


0, 


,42 



to 
CO 



»X10^® ergs s-1 in 3 - 80 keV 



Table 4. The best-fit PLCUTxCYAB2xGABS model parameters, determined by the flux-sorted spectra. 



flux 
level 



PLCUT 
Ecut 
(keV) 



-Bfold 
(keV) 



(keV) 



CYAB 



(keV) 



-Bo 2 
(keV) 



W2 
(keV) 



D2 



(keV) 



GABS 



(keV) 



Tga 



fdl 
fd6 



-0.34 
-0.26 



+0.09 
0.07 
+0.14 
0.07 



2.24 

1. 



+ 1.24 
1.24 
+ 1.95 
-0.69 



q + 0.73 
^-0.53 
g + 2-13 
"-0.80 



22.6 
25.3 



+0.6 
3.1 
+0.5 
1.0 



7.75 
8.44" 



+8.59 
1.03 
+ 1.22 
0.69 



1.03 
1.38 



+0.17 
0.70 
+0.17 
0.22 



+ 0.9 
^-0.9 
c + 1.0 
°-0.8 



+ 4.10 
'^-2.10 



7+10.61 
'3.71 
fi + 13.44 
°7.58 



27 24+° ™ 
29.55_Q 32 



3.20 
2.33 



+ 1.47 
-1.01 
+ 0.63 
0.76 



0.40" 
0.87 



+ 1.08 
-0.24 
+ 0.50 
0.39 



35.42 
14.50 



0.74 
0.61 



in 3 - 80 keV 
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Fig. 1. — Background-subtracted count rates and hardness ratios of X0331+53 obtained 
from the RXTE observations. Dashed hnes indicate the days of periastron passage (Zhang 
et al. 2005, and references therein), (a) The 2 — 10 keV ASM hghtcurve, with 1 day binning, 
(b) The 3-20 keV PCU2 hghtcurve. (c) The 20 - 80 keV HEXTE cluster B hghtcurve. (d) 
The hardness ratio, defined as the ratio between the HEXTE and PCU count rates. 
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Fig. 2. — Representative pulse-phase-averaged spectra of X0331+53, obtained on Dec 2a, 
Dec 24, Jan 20, and Feb 13b. (a) The background-subtracted and response-inclusive PCA 
(3 - 20 keV) and HEXTE (20 - 80 keV) spectra, (b) The PCA and HEXTE spectra in panel 
(a), normalized to the best-fit NPEX models which incorporates no absorption factors, (c) 
The data to NPEXxCYAB2 model ratios, (d) The data to NPEX xCYAB2x GABS model 
ratios. Detailed model parameters are given in Table [2J 



-27- 




5 10 20 50 100 

Energy (keV) 



Fig. 3. — (a) The same Dec 24 spectrum as shown in Figure [2|l, but presented in the 
deconvolved vF^, form. The filled circles represent the data, and the solid lines are the best- 
fit NPEXxCYAB2xGABS model. The dotted curve represent the best-fit model, from which 
the GABS absorption factor is removed. The dashed lines display the NPEX continuum, 
obtained by further resetting the CYAB2 factor to zero, (b) The best-fit model in its incident 
form. Meanings of the solid, dotted, and dashed curves are the same as in panel a. 
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Fig. 4. — Confidence contours between the continuum and tlie CRSF parameters, obtained 
from spectral fitting to the Dec 24 data. The solid, dashed, and dotted lines represent 
68%, 90%, and 99% confidence levels respectively. The crosses indicate 90%-confidence error 
ranges of the individual parameters given in Table 2. 
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Fig. 5. — The fundamental (panel a) and second harmonic (panel b) cyclotron resonance 
energies, shown against the 3 — 80 keV luminosity. Panel (c) shows the ratios between 
panels (b) and (a). The plotted data are extracted from the daily-averaged data using 
the NPEXxCYAB2xGABS model fitting. The error bars represent 90% confidence levels. 
Those data in which £'a2 is fixed at 2i?ai are omitted. The red and blue data represent the 
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Fig. 6.— The whole PCU 2 hghtcurve of X0331+53, plotted with 64 sec binnings. The 
horizontal dot lines represents boundaries of the flux sorting. The vertical dashed line divides 
the ascent and descent phases of the outburst. 
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they are shifted vertically by a factor of 0.5 between neighboring flux levels, (b) The spectra 
normahzed to the best-fit NPEXxCYAB2xGABS model. The data are shifted vertically, 
in the same way as in panel (a). 
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Fig. 8. — (a) Same as Figure [TJ but the data are derived from the outburst descent phase. 
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Fig. 9. — Same as Figure |5l but the results are derived using the flux-sorted spectra. The 
inset panel in (a) describes the luminosity dependence of the i?ga parameter. 
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Fig. 10. — The luminosity dependence of ^ai and ^a2- 
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Fig. 11. — (a) The flux-sorted fdl (red circles) and fd6 (blue triangles) spectra, normalized 
to the best-fit NPEX continuum determined individually using the 3-13 and 60-80 keV 
ranges only. The energy scale of the blue spectrum is scaled by a factor of 0.885, and it 
is compressed in the vertical direction by a factor of l+(l-y)x0.38 (with the ratio 1.0 kept 
as a pivot). The < 20 and > 20) keV ranges are covered by the PCA and HEXTE data, 
respectively, (b) The same as panel (a), but the PLCUT model is used as the continuum 
instead of NPEX. Again, the two continua are normalized to separate models. With respect 
to the red data points, the blue ones are scaled by a factor of 0.895 in the energy axis, and 
l+(l-y)x0.40 in the vertical direction. 
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Fig. 12. — The CRSF depths derived from the flux-sorted spectra, plotted against the 3 — 80 
keV source luminosity. 
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Fig. 13. — The estimated height of the CRSF forming regions against the 3 — 80 keV 
luminosity. 
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Fig. 14. — Comparison between X0331+53 and 4U 0115+63, in terms of the luminosity 
dependence of the fundamental resonance energy (panel a) and of the estimated cyclotron 
"photosphere" height. The CRSF are modeled with eqEl The results of 4U 0115+63 refer 
to Nakajima et al. (2006a). 
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Fig. 15.— CRSF depths D of X0331+53 plotted against the fractional width W/Ea. Results 
on 4U 0115+63 (INakajima et al. Il2006al ) are also shown for comparison. 



